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Abstract 

This  paper  presents  a  new  lattice-gas  method  for  molecular  dynamics 
modeling.  A  mean  field  treatment  is  given  and  is  applied  to  a  linear  stabil¬ 
ity  analysis.  Exact  numerical  simulations  of  the  solid-phase  crystallization 
is  presented,  as  is  a  finite-temperature  multiphase  liquid-gas  system.  The 
lattice-gas  method,  a  discrete  dynamical  method,  is  therefore  capable  of 
representing  a  variety  of  collective  phenomena  in  multiple  regimes  from 
the  hydrodynamic  scale  down  to  a  molecular  dynamics  scale. 
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1  INTRODUCTION 

Presented  in  the  paper  is  a  theory  of  lattice-gas  dynamics  that  includes  inter¬ 
particle  potentials.  The  microscopic  lattice-gas  dynamics  is  a  highly  discrete 
form  of  traditional  molecular  dynamics.  All  the  usual  dynamical  quantities  ap¬ 
pearing  in  the  traditional  theory  are  discrete  in  the  case  of  a  lattice-gas.  It  is 
well  known  that  in  lattice-gases  these  discrete  quantities  include  space,  time, 
and  momentum.  Here  the  notion  of  a  discrete  field  is  introduced  and  an  analyt¬ 
ical  theory  is  presented  to  describe  emergent  macroscopic  dynamics.  The  use  of 
the  field  concept  is  shown  to  be  quite  useful.  In  particular,  the  field  concept  is 
useful  in  describing  our  novel  lattice-gas  with  multiple  long-range  interactions 
with  different  ranges  and  polarity.  This  new  lattice-gas  possesses  a  liquid-solid 
transition  and  can  be  used  as  a  new  general  method  of  simulating  molecular 
dynamics.  The  theoretical  possibilities  for  such  a  lattice-gas  opens  the  subject 
of  exactly  computable  modeling  to  the  areas  of  dynamical  solid-state  systems. 

It  is  known  that  interparticle  potentials  can  be  modeled  by  including  a  sin¬ 
gle  anisotropic  fixed-range  interaction  in  the  lattice-gas  dynamics  for  discrete 
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momentum  exchange  between  particles.  The  simplest  theoretical  model  of  this 
kind  is  the  Kadanoff-Swift-Ising  model  [8] .  An  attractive  fixed-range  interaction 
was  used  in  a  lattice-gas  automaton  by  Appert  and  Zaleski  [2]  in  1990  to  model 
a  nonthermal  liquid-gas  phase  transition.  The  use  of  attractive  and  repulsive 
fixed-range  interactions  of  this  sort  extended  the  lattice-gas  dynamics  to  a  fi¬ 
nite  temperature  liquid-gas  transition  where  a  complete  pressure,  density,  and 
temperature  equation  of  state  is  modeled,  and  the  complete  liquid-gas  coexis¬ 
tence  curve  is  analytically  predicted  through  a  Maxwell  construction  [15].  Our 
finite  temperature  liquid-gas  lattice-gas  is  presented  in  the  paper  for  pedagogi¬ 
cal  reasons  as  well  as  to  validate  the  theoretical  method  presented.  Lattice-gas 
crystallization  is  introduced  as  a  direct  generalization  of  the  finite  temperature 
liquid-gas  lattice-gas  model. 

This  paper  is  organized  in  basically  two  parts.  The  first  part  of  the  pa¬ 
per  up  to  and  including  section  §4  is  well  known  by  the  lattice-gas  community 
and  is  given  here  as  review  material  for  the  subject  of  lattice-gases  with  purely 
local  collisions.  The  rest  of  the  paper  presents  new  results,  using  a  new  theo¬ 
rem  referred  to  here  as  the  lattice  multiple  theorem  presented  in  appendix  12. 
This  theorem  is  useful  for  determining  the  linear  response  of  a  lattice-gas  with 
long-range  interactions.  Appendix  12  describes  in  some  detail  an  explicit  nu¬ 
merical  method  for  implementing  the  the  simplest  of  long-range  interactions: 
the  bounce-back  and  clockwise  orbits. 


2  Lattice-Gas  Automaton:  An  Exactly  Computable 
Dynamical  System 

A  boolean  formulation  of  an  exactly  computable  dynamical  system,  known  as 
a  lattice-gas,  may  be  stated  in  a  way  that  is  consistent  with  the  Boltzmann 
equation  for  kinetic  transport.  In  essence  the  lattice-gas  dynamics  are  a  sim¬ 
plified  form  of  molecular  transport  as  we  restrict  ourselves  to  a  discrete  cellular 
phase  space.  The  macroscopic  equations,  in  particular  the  continuity  equation 
and  the  Navier-Stokes  equation,  are  obtained  by  coarse-graining  over  a  discrete 
microdynamical  transport  equation  for  number  boolean  variables.  The  scheme 
employs  the  finite-point  group  symmetry  of  a  crystallographic  spatial  lattice.  It 
is  somewhat  inevitable  that  to  obtain  an  exactly  computable  representation  of 
fluid  dynamics  one  must  perform  a  statistical  treatment  over  discrete  number 
variables. 

Before  introducing  the  basic  lattice-gas  microdynamical  transport  equation, 
let  us  give  some  notational  conventions.  We  consider  a  spatial  lattice  with  N 
total  sites.  The  fundamental  unit  of  length  is  the  size  of  a  lattice  cell,  I,  and 
the  fundamental  unit  of  time,  r,  is  the  time  it  takes  for  a  speed-one  particle  to 
go  from  one  lattice  site  to  a  nearest  neighboring  site.  Particles,  with  unit  mass 
m,  propagate  on  the  lattice.  The  unit  lattice  propagation  speed  is  denoted 
by  c  =  K  Particles  occupy  this  discrete  space  and  can  have  only  a  finite  B 
number  of  possible  momenta.  The  lattice  vectors  are  denoted  by  Cai  where 
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a  =  For  example,  for  a  single-speed  gas  on  a  triangular  lattice, 

a  =1,2,. ..,6.  A  particle’s  state  is  completely  specified  at  some  time,  t,  by 
specifying  its  position  on  the  lattice,  Xi,  and  its  momentum,  pi  =  mcCai,  at  that 
position.  The  particles  obey  Pauli  exclusion  since  only  one  particle  can  occupy 
a  single  momentum  state  at  a  time.  The  total  number  of  configurations  per  site 
is  2^.  The  total  number  of  possible  single  particle  momentum  states  available 
in  the  system  is  A’totai  =  BN.  With  P  particles  in  the  system,  we  denote  the 
filling  fraction  hy  d  =  ■ 

The  number  variable,  denoted  by  na{x,  t),  takes  the  value  of  one  if  a  particle 
exists  at  site  x  at  time  t  in  momentum  state  mcBa,  and  takes  the  value  of  zero 
otherwise.  The  evolution  of  the  lattice-gas  can  then  be  written  in  terms  of  as 
a  two-part  process:  a  collision  and  streaming  part.  The  collision  part  reorders 
the  particles  locally  at  each  site 

n'„(x,t)  =  na(x,t) -I- na(n(x,t)),  (1) 

where  fla  represents  the  collision  operator  and  in  general  depends  on  all  the 
particles,  ft  at  the  site.  So  as  a  short-hand  we  suppress  the  index  on  the  occu¬ 
pation  variable  when  it  is  an  argument  of  fla(ji{x,t))  to  represent  this  general 
dependence.  In  the  streaming  part  of  the  evolution  the  particle  at  position  x 
“hops”  to  its  neighboring  site  at  x  -|-  leg^  and  then  time  is  incremented  by  t 

n'„(x  -h  le^,  t  +  T)  =  na(x,  t)  +  fIo(n(x,  t)).  (2) 

Equation  (2)  is  the  lattice-gas  microdynamical  transport  equation  of  motion. 
The  collision  operator  can  only  permute  the  particles  locally  on  the  site  since 
we  wish  the  local  particle  number  to  be  conserved  before  and  after  the  collision. 

We  construct  an  n-th  rank  tensor  composed  of  a  product  of  lattice  vectors 
[13] 

=^(ea)n---(ea)i„,  (3) 

a 

where  a  =  1, . . .  ,B.  All  odd  rank  E  vanish.  We  wish  to  express  in  terms 

of  Kronecker  deltas,  6ij  =  1  for  i  =  j  and  zero  otherwise.  We  can  turn  this 
problem  of  expressing  the  A-tensors  in  terms  of  products  of  Kronecker  deltas 
into  a  problem  of  combinatoric  counting.  We  use  the  following  tensors 

A?,  =  %  (4) 

‘^ijki  ~  dijSki  +  dikSji  +  SiiSkj  (5) 

and  so  forth.  Then  we  know  that  if  E  is  isotropic  it  must  be  proportional  to  A 

E{2n)  ^  ^(2«) 

and  that  the  constant  of  proportionality  may  be  obtained  by  counting  the  num¬ 
ber  of  ways  we  could  write  a  term  comprising  a  product  of  n  Knonecker  deltas. 
Consider  for  example  the  case  n  =  2.  Since  the  Knonecker  delta  is  symmetric  in 
its  indices,  the  following  four  products  are  identical:  SijSki  =  Sijdik  =  SjiSki  = 
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^jiSik-  The  degeneracy  is  2^.  Furthermore,  the  order  of  the  Kronecker  deltas 
also  doesn’t  matter  since  they  commute;  that  is,  Sijdki  =  SkiSij.  The  degen¬ 
eracy  is  2!.  For  the  case  where  n  is  arbitrary,  there  are  2”  identical  ways  of 
writing  the  product  of  n  Kronecker  deltas.  For  each  choice  of  indices,  there  are 
an  additional  n!  number  of  ways  of  ordering  the  products.  Therefore,  the  total 
number  of  degeneracies  equals  2"n!  =  (2n)!!.  The  total  number  of  permutations 
for  2n  indices  equals  (2n)!.  So  from  this  counting  procedure  we  know  that 
consists  of  a  sum  of  ^yyr  =  (2n  —  1)!!  terms. 

In  general,  the  lattice  tensors  are 


^2n+l  ^  Q 

Jp‘2n  _ ^ ^  2n 

D{D  +  2)---{D  +  2n-2) 


(7) 

(8) 


3  Coarse-Grained  Dynamics 

To  theoretically  analysize  the  lattice-gas  dynamics,  it  is  convenient  to  work  in 
the  Boltzmann  limit  where  a  field  point  is  obtained  by  an  ensemble  average 
over  the  number  variables.  That  is,  we  may  define  a  single  particle  distribution 
function,  fa  =  (u-a),  resulting  from  an  ensemble  of  initial  conditions  and  the 
neglect  of  correlations,  with  the  averages  taken  over  the  ensemble. 

It  is  essential  to  determine  the  macroscopic  limit  of  the  microdynamical 
transport  equation  (2)  and  to  see  how  it  leads  to  non-compressible  viscous 
Navier-Stokes  hydrodynamics — for  a  lengthier  treatment  of  this  see  Frisch  et 
al.  [6]. 

Using  the  Boltzmann  molecular  chaos  assumption  the  averaged  collision  op¬ 
erator  simplifies  to  {fla{fi))  =  Ua((n)),  and  by  coarse-graining  and  Taylor  ex¬ 
panding  (2)  we  obtain  the  lattice  Boltzmann  equation 

dtfa  +  CCaidifa  =  ^a{f)-  (9) 

We  write  the  particle  number  density,  momentum  density,  and  moment  flux 
density  in  terms  of  the  single-particle  distribution  function  as  follows 

fa  =  P 
a 

mc^ea^fa  =  pVi 
a 

‘ITIC  ^  ^  GaiCqj  fa  —  Tlij 
a 

Now  following  Landau  and  Lifshitz  [9]  we  know  that  in  standard  form  we 
must  be  able  to  write  the  momentum  flux  density  tensor  as  follows 

mc^  ^  CaiCajfa  =  pS^j  +  pViVj  -  (13) 


(10) 

(11) 

(12) 
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where  in  (13)  the  first  two  terms  represent  the  ideal  part  of  the  momentum  flux 
density  tensor  and  =  r]{diVj  —  djVi)  is  the  viscous  stress  tensor.  Alternatively 
the  momentum  flux  density  tensor  may  be  written 

Hy  =  mc^  ^  eaiCajfa  =  -(Tij  +  pViVj,  (14) 

a 

where  cr^  is  the  pressure  stress  tensor 

CTy  =  -pSij  +  7]  (diVj  +  diVj) .  (15) 


The  general  form  of  the  single  particle  distribution  function,  appropriate  for 
single  speed  lattice-gases,  is  a  Fermi-Dirac  distribution.  Fundamentally,  this 
arises  because  the  individual  digital  bits  used  to  represent  particles  satisfy  a 
Pauli-exclusion  principle.  Therefore,  the  distribution  must  be  written  as  a  func¬ 
tion  of  the  sum  of  scalar  collision  invariants,  a  +  iScaiVi,  implying  the  following 
form  ^ 

1  _|_  ga+0eaiVi  ' 

Taylor  expanding  (16)  about  u  =  0  to  fourth  order  in  the  velocity  and  equating 
the  zeroth,  first,  and  second  moments  of  fa  to  (10),  (11),  and  (12)  respectively, 
the  parameters  a  and  /3  are  determined.  The  inviscid  part  of  the  lattice-gas 
distribution  function  becomes 


/  £  eq\ideal 

)lga 


n 

B 


nD  nD{D  +  2)  ^  ^ 

^ai^i  9  2,C^ B 


ti{D  +  2)  2 


(17) 


where 


D  1  -  2d 
^  ^  D  +  2  1-d' 


(18) 


That  is,  using  p  =  mn  for  the  density  and  Cs  =  for  the  sound  speed,  the 
moments  of  lattice-gas  distribution  are 


(  f  <=<lYd.eal 
^  2^\Ja  )lGA 

=  p 

(19) 

rnc  2_^  fiai\Ja  >  LG  A 
a 

=  pVi 

(20) 

rnc  2_^taitaj{Ja  )lGA 

=  pc2(l-g— )% 

+gpViVj. 

(21) 

The  lattice-gas  automaton  almost  produces  the  correct  form  for  the  momentum 
flux  density  tensor,  except  that  11^  appears  to  have  a  spurious  dependence  on 
the  square  of  the  velocity  field,  (1  —  g^)  with  a  factor  g  arising  as  an  artifact 
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of  the  discreteness  of  the  number  variables.  Working  directly  in  the  Boltzmann 
limit  and  using  only  symmetry  arguments,  it  is  possible  to  fix  this  problem. 

The  macroscopic  equations  of  motion  are  then  determined  from  mass  con¬ 
servation  (continuity  equation)  and  momentum  conservation  (Euler’s  equation) 

dtp  +  d^{f)Vi)  =  Q  (22) 

and 

dt{pvi)  +  djiiij  =  0.  (23) 

Substituting  (21)  into  Euler’s  equation  (23),  gives  us  the  Navier-Stokes  equation 
for  a  viscous  fluid 

p  {dm  +  gvjdjVi)  =  -dip  +  pd'^v^,  (24) 

given  a  non-divergent  flow  {diVi  =  0)  appropriate  to  the  incompressible  fluid 
limit  and  where  the  pressure  is 

p  =  pci  (^1  +  .  (25) 

A  general  expression  for  the  shear  viscosity,  p,  for  a  single-speed  lattice-gas  has 
been  derived  by  Henon  [7]. 

In  any  lattice-gas  simulation,  one  typically  obtains  a  realization  of  the  macro¬ 
scopic  dynmical  variables  by  block  averaging  in  both  space  and  time  over  the 
microscopic  variables.  In  this  way,  for  example,  a  momentum  map  can  be  pro¬ 
duced  so  that  the  dynamic  evolution  the  the  fluid  can  be  monitored.  The  size 
of  the  coarse  grain  block  affects  the  resolution  with  which  one  can  observe  the 
system  but  of  course  does  not  at  all  affect  the  underlying  dynamics.  If  too 
small  of  a  coarse  grain  block  size  is  used,  more  fluctuations  in  the  macroscopic 
variables  occurs. 


4  Lattice  BGK  Equation 


We  wish  to  consider  a  dynamical  transport  equation  for  the  particle  distribution 
function  given  in  the  previous  section.  We  have  a  lattice  Boltzmann  gas  defined 
on  a  discrete  spatial  lattice.  Restricting  ourselves  to  a  single  speed  lattice-gas 
system,  the  lattice  BGK  equation  is 

^  -k  cea^d^fa  =  -^{fa-  /a'"'')-  (26) 

This  equation  was  introduced  in  1954  by  D.  Bhatnager,  E.  Gross,  and  M.  Krook 
[3,  10].  A  way  to  obtain  (26)  was  introduced  by  Shiyi  Ghen  et  al.  [5]  by  ex¬ 
panding  the  lattice  Boltzmann  collision  term  to  first  order  about  the  equilibrium 
distribution  and  assuming  it  diagonal. 

It  is  possible  to  fix  the  anomaly  in  the  fluid  pressure  that  occurs  in  the 
lattice-gas  automaton.  Ghen  et  al.  [4]  have  introduced  a  pressure-corrected 
equilibrium  distribution  to  have  the  following  Ghapman-Enskog  expansion 


(  f  eq\ideal  _  ^ 

yJa.  )bgk  —  ^ 


nD 


^ai^i  “k 


nD{D  +  2) 
2c2R 


nD 

2c^B^ 


(27) 
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which  satisfies 


=  p(i-w) 

(28) 

^  ggz  \  Ja  )  PC 

II 

(29) 

p  -p  (  f 

=  pclS^j  +  pViVj. 

(30) 

a 
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Here  the  definition  of  the  density  is  modified  by  the  1  —  ^  factor. 

5  Density  Dependent  Pressure  in  the  Boltzmann 
Limit 

Here  we  exploit  the  analytical  facility  of  the  lattice  Boltzmann  approach  and 
show  that  the  addition  of  a  convective-gradient  term  in  the  lattice  Boltzmann 
equation  allows  one  to  model  a  hydrodynamic  gaseous  flow  governed  by  a  general 
equation  of  state  [14].  The  pressure  may  have  a  nonlinear  dependence  on  the 
local  density.  It  is  possible  to  generalize  this  to  a  multi-speed  lattice-gas  or  a  to 
single-speed  lattice-gas  coupled  to  a  heat  bath  so  that  the  pressure  dependence 
includes  the  local  temperature  as  well. 

The  equation  of  state  for  the  isothermal  gas  is 

P  =  clp.  (31) 

We  now  wish  to  consider  how  we  may  alter  the  lattice-Boltzmann  equation  to 
allow  for  a  more  general  equation  of  state.  Let  us  add  an  additional  term, 
Pa{x,  X  +  rca),  to  the  R.H.S.  of  (9) 

dtfa{x,t)  +  ceaidifa{x,t)  =  - [fla (a;,  t) -f 

T 

Pa{x,x  +  rea,t)].  (32) 

Pa  depends  on  the  local  configuration  of  the  system  at  position  x  as  well  as  on 
the  local  configuration  at  a  remote  position  x  +  rCa-  We  assume  that  the  values 
oiif)  dX  X  and  x  +  rCa  are  independent  and  that  therefore  Pa  can  be  factorized  ^ 

Pa{x,x  +  rea,t)  = '>p{x)2p{x  +  rCa).  (33) 

^The  effect  of  long-range  interactions  on  fa{x)  will  actually  depend  on  local  configurations 
at  fc  +  rca  and  at  x  —  rca.  So  we  could  write  the  full  form  of  the  long-range  part  of  the  collision 
operator  as  ^  rCa)  —  —  'rea)V’(^)]i  where  th^  factor  of  ^  must  be  included  to 

avoid  double  counting  when  doing  any  directional  sums,  since  ^  here  does  not  have 

any  directional  dependence.  According  to  Theorem  1  in  the  appendix  12,  upon  expanding 
±  rea),  both  terms  would  add  to  remove  the  ^  factor,  so  using  Pa  =  'rp(x)'\l){x  +  rCa) 
in  the  present  calculation  ultimately  gives  the  same  result.  In  appendix  12  where  we  give  the 
microscopic  long-range  part  of  the  lattice-gas  collision  operator,  the  microscopic  -0  does  have 
directional  dependence  so  we  use  the  full  form  there. 
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In  a  single  speed  lattice-gas  model  as  we  have  been  considering,  Pa  is  a  function 
of  the  local  density.  In  a  single-speed  model  coupled  to  a  heat  bath,  Pa  may 
depend  on  the  local  temperature  as  well  [15]. 

We  wish  to  constrain  the  form  of  Pa  so  as  not  to  violate  continuity.  We 


require 

o' 

II 

(34) 

and  when  9^/^  =  0  , 

E  =  0. 

(35) 

a 


Constraint  (35)  is  required  only  under  uniform  filling  conditions;  i.e.  for  general 
situations  e-aiPa  is  non-zero.  In  the  uniform  flow  limit  the  lattice-Boltzmann 
equation  reduces  to 

dtfa{x,  t)  =  -  [f^a(a;,  t)  +  P{x,  x  +  rca,  t)] ,  (36) 

T 

where  we  have  taken  the  directional  dependence  of  the  long  range  collisional 
term  to  occur  only  in  its  argument,  Pa{x  +  rca)  — >■  P{x  +  rCa)-  We  assume 
the  probability  of  a  long  range  collision  depends  only  the  density  at  the  spatial 
location  of  a  momentum  transfer  event  and  not  on  the  direction  of  the  momen¬ 
tum  transfer.  That  is,  we  require  the  interaction  distance  to  be  of  sufficiently 
long  range  that  the  approximation  of  local  isotropy  in  the  particle  distribution 
is  valid.  Summing  over  all  lattice  directions  and  using  constraint  (34)  we  have 
maintained  the  collision  property  that 

=  =  (37) 

a  a 

Thus,  for  arbitrary  flows,  summing  the  lattice-Boltzmann  equation  (32)  over  all 
directions  preserves  continuity 

=  0  (38) 

=  0,  (39) 


dt'^faP  c'^Caidifa 

a  a 

— ^  dtp+di{pv,) 


where  we  have  used  (37). 

Multiplying  the  lattice  Boltzmann  equation  by  Cai  and  then  summing  over 
directions  gives 


^  (  ^aifa{x^  t)  cdj  ^  f a{,X ■,  t)  — 


-'Ipix,  t)  E  eazipix  +  rCa,  t). 


(40) 


Using  Theorem  1  in  appendix  12  we  can  expand  the  R.H.S. 


(I)  "  i) 

(41) 


Therefore,  we  again  arrive  at  Euler’s  equation 


dt{pu,)+dj{ll,j)  =  0  (42) 

but  with  an  augmented  momentum  flux  density  tensor 

=  mc^^eaieajfa-mc^B  /  ^Xk  ^ix,t)dk{rd)~^ I D{rd)'>p{x,t). 

(43) 

Since  the  additional  term  in  the  momentum  flux  density  tensor  is  diagonal,  it 
can  only  impart  an  effective  density  dependent  pressure. 

Defining  a  configurational  potential  energy  as 

E(x)  = 

X  J  dxk  ip{x,t)dk{rd)~^ lD{rd)tp{x,t)  (44) 

then  Euler’s  equation  (42)  gives  us  the  viscous  Navier-Stokes  equation  for  non¬ 
ideal  fluids 


dt{pv,)  +  djipviVj)  =  -d,  {cIp  +  V (p))  +  pvd'^v^.  (45) 

Therefore,  we  have  arrived  at  a  general  equation  of  state  defined  by  the  potential 
energy  function  V{p)  where  there  is  an  inter-particle  force  Ei(x)  =  —diV{p{'x)). 
The  form  of  the  density  dependent  pressure  directly  follows 

P{p)  =  clp+V{p).  (46) 

With  this  methodology,  we  can  model  a  system  with  a  general  equation  of  state 
with  completely  local  dynamics  described  by  the  generalization  of  (26) 

fa{x  +  le,t  +  T)  =  faix,t)  -  ^{fa{x,t)  -  fa''\x,t)) 

+'ip{x,t)'ip{x  +  rea,t).  (47) 

In  the  Boltzmann  limit,  the  analysis  itself  does  not  indicate  the  form  of  V  in 
(46)  (or  more  to  the  point,  does  not  indicate  the  form  of  ip),  but  does  show  it  is 
possible  to  have  a  lattice-gas  that  has  Navier-Stokes  dynamics  as  its  macroscopic 
limit  with  a  density  dependent  pressure  (45).  This  is  the  motivation  needed  to 
develop  a  more  complete  microscopic  description.  With  a  lattice-gas  automaton 
microscopic  description,  the  interparticle  force  —diV{p)  may  be  caused  by  long- 
range  momentum  exchange  between  two  particles.  Calculating  the  probability 
of  such  momentum  exchange  events  should  provide  a  way  to  determine  tp. 

Note  that  in  the  mesoscopic  regime  in  which  the  Boltzmann  equation  is 
applicable,  the  lowest  order  expression  for  V  is  proportional  to  That  is 

V{x)  =  J  dxk'ip{x)dk'ip{x)  (48) 
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or 


(49) 


mc^B  /r\  2/  ^ 

A  similar  calculation  has  also  been  done  by  Shan  and  Chen  [12]  who  have  verified 
their  analysis  by  comparing  with  data  taken  from  lattice  BGK  simulations. 
They  have  also  presented  exact  calculations  for  the  liquid-gas  interface  profile 
and  surface  tension.  In  the  following  section  §6  we  take  another  view  point 
and  write  an  alternate  expression  for  the  potential  energy,  but  one  that  is  also 
proportional  to  ijP'  .  This  alternate  view  of  the  potential  energy  will  help  towards 
developing  the  lattice-gas  automaton  microscopic  description. 

6  Interaction  Energy 

We  introduce  a  potential  energy  due  to  non-local  2-body  interactions 

'  =  ^  X!  X!  fo.{x))fb{x)Vabmn{l-  fm{x))fn{x),  (50) 

{xx')  {abmn) 

where  Vabmn  =  V ^abmn  and  Aabmn  IS  either  ±1  or  vanishes  for  any  set  {abmn} 
that  violate  mass,  momentum,  or  energy  conservation.  H'  accounts  for  the 
potential  energy  between  particles  in  coming  along  lattice  directions  b  and  n 
and  outgoing  along  a  and  m  and  is  therefore  restricted  to  2-body  interactions. 

We  now  try  to  justify  the  form  of  H'  and  in  so  doing  develop  a  microscopic 
description  of  the  lattice-gas  with  long  range  interactions.  We  require 

Pi  =  -d,H'.  (51) 

H'  is  thought  of  as  the  configurational  potential  energy  due  to  momentum  trans¬ 
fers  between  two  locations.  The  momentum  exchange  per  unit  time  between  two 
points  X  and  x'  in  the  fluid  is 

5pi  =  —  mvf,  (52) 

where  the  incoming  and  outgoing  velocity  states  are  quantized:  vf  =  ceu  and 
yout  _  probability  '4’{x)  of  there  being  a  local  momentum  change  at 

some  point  x  depends  independently  on  the  probability  fb{x)  that  there  is  a 
particle  in  velocity  state  ceu  and  the  probability  1  —  fa{x)  there  isn’t  a  particle 
in  velocity  state  ccai-  So  in  this  factorized  approximation  that  neglects  particle- 
particle  correlations,  we  write 

Tpix)  =  (1  -  fa{x))fb{x).  (53) 

As  a  long  range  momentum  exchange  event  involves  two  sites,  x  and  x' ,  we  can 
define  the  vector  =  rcai  =  Xi  —  x[  and  the  therefore  parallel  and  perpendicular 
components  of  the  local  momentum  exchange  are 

5p\\  =  5p  •  f 

6p±  =  I  (5p  X  f  I  . 
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(54) 

(55) 


The  two  components  of  the  force  mediated  by  the  long  range  momentum  ex¬ 
change  could  be  interpreted  as  created  by  two  separate  fields 


(5E 

(5B 


Sp\\{r)c^ 


6p±{r)c 

I 


5p  X  f , 


(56) 

(57) 


where  6p\\  and  Sp±  are  written  as  a  function  of  r  since  any  kind  of  functional 
dependence  is  allowed  provided  enough  detail  is  specified  for  the  automaton 
interaction  rules.  We  have  explicitly  written  the  forms  of  the  parallel  and  per¬ 
pendicular  components  of  a  lattice-gas  force  field  to  stress  an  analogy  with  the 
classical  theory  of  electromagnetism  where  the  electric  and  magnetic  fields  are 
expressed,  for  a  differential  element  of  charge  and  current  element  respectively, 
by  the  well  known  laws  of  Coulomb  and  Biot-Savart 


dE 

dB 


SQ  . 


47r7-2 

I 

dTrr^ 


r 

dl  X  f . 


(58) 

(59) 


Using  (53)  we  can  write  the  total  momentum  change  as  a  field  itself 


Spi(^X^  X  )  —  TTIC  ^  ^  (Cai  6nz)(l  ^a(^))^fc(^)-^a6mn(l  fmi,^  )■ 

{ahmn) 

(60) 

Note  that  the  momentum  change  is  zero  for  a  fluid  with  uniform  density  due 
to  the  symmetry  of  the  lattice.  For  central  body  interparticle  momentum  ex¬ 
changes,  in  the  Boltzmann  limit  we  can  then  approximate  the  configurational 
potential  energy  as 

E(a;,a:')  = (y)  ^  r-(ea-e6-kem-e„)(l-/a)/&Aabm„(l-/m)/n, 

{abmn) 

(61) 

where  r  is  the  range  of  the  interaction.  For  a  system  locally  isotropic  in  its 
particle  distributions,  letting  i/)(a;)  =  (1  —  f{x))f{x),  this  may  be  simplified  to 

V{x,x')  =  —-mc^  (y  )  E  E  r-(ea-e{,-ke 

m  )Kbmni^{x)tp{x'),  (62) 

{abmn)  {abmn) 


which  is  suitable  for  a  bulk  description  of  the  fluid.  Now  the  form  of  H'  follows 
if  we  sum  over  all  pairs  {xx')  and  define 

^abmn  —  XflC  (y)  ^  '  (®<i  -f-  ^n^^abmn^  (63) 

SO  that 

H'  =  '^  V{x,x')  =  Y  '^i.x)Vabmn'ip{x').  (64) 

{xx')  {xx')  {abmn) 
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7  E  Field  Construction 


It  is  possible  to  define  a  field  that  exists  in  a  lattice-gas  that  has  long-range 
momenta  exchanges  occurring.  The  notion  is  to  consider  each  lattice-gas  particle 
as  having  a  delta  function  type  field  that  exists  only  at  certain  fixed  ranges  and 
certain  fixed  angles.  Therefore,  the  lattice-gas  particle  has  a  highly  anisotropic 
field.  However,  in  the  coarse-grained  limit  obtained  by  averaging  over  many 
particles,  a  valid  description  of  a  continuous  field  emerges.  Of  coarse,  if  there 
are  no  gradients  in  the  coarse-grained  density  of  the  system,  the  field  must 
necessarily  vanish. 

The  field  at  position  x  due  to  a  particle  at  position  x'  along  the  lattice 
direction  a  is  a  delta  function 

eaiix]  x')  =  ^  a{a)S{x  -  x'  -  r„ea)eai  (65) 

(7 

That  is,  the  discrete  field  must  be  directed  along  and  must  be  a  distance 
from  the  source,  ov  x  =  x'  +  r^Ca-  The  total  field  is  obtained  by  considering  all 
the  possibilities  where  a  particle  could  contribute.  The  sum  over  a  is  necessary 
to  account  for  multiple  interaction  ranges,  where  Q;(cr)  denotes  the  strength  of 
the  interactions  at  range  r^.  Thus  to  obtain  the  total  field  we  must  sum  over 
all  directions  and  integrate  over  all  positions 

Ei{x)  ='^  f  dx'^|){x')ea^{x;x')  (66) 

a 

or 

TTIC^ 

Ejjx)  =  —i—y^a{a)y^eai^{x  -  rgeg).  (67) 

(7  a 

Using  Theorem  1  we  can  evaluate  the  directional  sum  and  express  the  field  as 
a  gradient  of  a  scalar  quantity 


Et{x)  =  -di 


(I)'  {r^d) 

<7 


?  lD{r^d)'tp{x)  , 


(68) 


where  is  defined  as 

=  -a(a)  ■  (69) 

Note  that  /x  >  0  for  attractive  interactions  and  /i  <  0  for  repulsive  interactions. 
Since,  Ei  =  di(j),  the  field’s  scalar  potential  is 


(j)  =  (r^d)  ?  I D{r^d)'ijj{x) 


(70) 


=  —mc^B^^i 


— U  H - 7 - rr'^d'^'tp  ■ 

2D{D  + 2)  ^ 


(71) 
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The  field  Ei  and  the  event  probability  ip  =  (i(l  —  d)  appear  in  the  Navier-Stokes 
equations  as  follows 


dtipvi)  +  dj{gpViVj) 


-cldi 


P 


+  ipEi  +  pvd'^Vi 


(72) 


and  according  to  (43),  the  term  ipE^  when  due  to  central-body  interactions  can 
modify  only  the  pressure. 


8  Stability  Analysis 


In  this  section  we  consider  the  linear  response  of  a  lattice-gas  with  long  range 
central-body  interactions.  The  macroscopic  equations  of  motion  are  (39),  (72), 
and  (68)  respectively 


^tp  +  ^i{pv^)  =  0 

dtipVi)  +  dj{gpviVj)  =  -cld. 


p[  1  +  ff- 


-I-  ipEi  +  pvd'^Vi 


Ei{x)  =  -di 


mc^B  {r^d)  ?  lD{r„d)ip{x) 


We  treat  the  effect  of  the  field  Ei  as  a  perturbation  on  a  resting  equilibrium 
state  where  p  is  uniform  and  constant  and  v  =  0.  Then  an  e-expansion  of  the 
dynamical  variables  is 


Vi  =  eui 
p  =  Po  +  £Q 
Ip  =  ipo  +  ep- 


Using  Ip  =  {1  —  d)d  we  have 


P  = 


1  -  2do 
mB 


where  do  =  Consequently,  the  linear  response  equations  are 
dtP  -f  Po^illi  —  0 

PodtUi  =  -cldip  -  ipodi  mc^B'^pa-  ""  irod)~^lD{rod)p 


(73) 

(74) 

(75) 

(76) 

(77) 

PoVod'^Ui 

(78) 


Then  applying  dt  to  the  continuity  equation  and  di  to  the  Navier-Stokes  equation 
allows  us  to  eliminate  Ui  and  to  obtain  the  following  second-order  equation  in  g 


dfP  =  cld'^g  +  ipo{l  -  2do)c^d‘^ 


(rod)  ’^lD{rod)g 


Vod'^dtp. 

(79) 
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In  an  inviscid  fluid  {v  =  0)  with  no  interparticle  potentials  {Ei  =  0)  g  would 
satisfy  the  wave  equation 


(80) 


Given  a  non-zero  perturbation,  g  can  be  Fourier  expanded 

g  =  j  dujdh  (81) 

and  we  can  replace  g  with  g  by  taking  dt  — >■  —ito  and  di  — >■  iki.  (79)  becomes 


—oj'^g  =  —k^clg—il^o{^  —  ‘^do)k^c^ 


(r^k)  ’^JD{rak)g 


-ioji'ok'^g, 
(82) 

where  we  have  made  use  of  the  identity  that  relates  the  hyperbolic  Bessel  func¬ 
tion  with  imaginary  argument  to  the  ordinary  Bessel  function 


(iz)  =  z  ''J„{z). 

Dividing  out  g  gives  a  quadratic  equation  for  uj 


(83) 


UJ  \  Vo 


—  +1- 


-\-k^ 


1-k 


V'o(l  -  2fio)7?^Ma  ""  (r^fc)  ^Joir^k) 


The  dispersion  relation  for  uj{k)  is  then 


(84) 


tpoi^  -  2do)D'^fi„  (yy  (rak)  ’^jD{r„k) 


2c, 


(85) 


In  the  long  wavelength  limit,  (85)  reduces  to  linear  sound  speed  dispersion 

UJ  =  Cak.  (86) 

In  the  absence  of  long  range  interactions,  (85)  reduces  to  the  dispersion  relation 
for  ideal,  incompressible,  viscous  flow 


UJ  =  Cgk 


v'^k'^ 

4c2 


(87) 


By  choosing  different  ranges  and  strengths  of  the  momentum  exchanges  we 
can  adjust  the  dispersion  relation  (85)  and  produce  a  variety  of  interesting  dy¬ 
namical  behavior.  Therefore,  in  the  macroscopic  limit,  what  principally  defines 
the  linear  response  of  a  lattice-gas  model  with  long  range  interactions  is  the  set 
of  constants  a(cr)  and  To-. 
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9  Finite  Temperature  Liquid-Gas  Model 

A  simple  two-dimensional  example  is  a  lattice-gas  on  a  triangular  lattice.  The 
macroscopic  equations  of  motion  for  the  lattice-gas  automaton  are 


dtp+di{pvi)  =  0 

dt{pvi)  +  dj{gpv^Vj)  =  -diP  +  pvd^Vi 


dip  =  cldi 


P[  1  +  5- 


-  i’Ei. 


The  field  to  lowest  order  in  the  interaction  range  is 

mc^B  /r' 


E^{x)  =  -- 

This  implies  that  the  local  force  is 


f/'Ai  =  -di 


D 


a 


(j) 


^-mc,B  I  T 


(7)*= 


Since  ip  =  d{l  —  d),  the  pressure  is  a  simple  polynomial 

a{h)r 


p{d,  h)  =  mc^B 


^  1  +  5-y  + 


21 


d\l  -  dy 


(88) 

(89) 

(90) 

(91) 

(92) 

(93) 


where  we  have  written  the  pressure  depending  on  the  density  and  a  temperature 
control  parameter,  h,  that  modifies  the  strength  of  the  interactions,  a  =  a{h). 
This  will  be  discussed  in  more  detail  below.  It  is  possible  to  determine  the 
coexistence  curve  for  such  a  finite  temperature  liquid-gas  system.  We  begin  by 
defining  the  free  energy,  G,  as 


G{d,h)=  f 

J  do 


dn  dp{n,  h) 
dn 


(94) 


Using  (93)  for  a  fluid  at  rest  we  can  carry  out  the  integral  to  obtain 


G{d,  h)  =  mclB 


x{h)i 


4d3 


4d„ 


—2  d  -|-  3  - ^ — h  2  do  —  3  do  H - ^ — J  -I-  log 

(95) 

A  Maxwell  construction  can  be  performed  by  making  a  parametric  plot  of  the 
free  energy  versus  the  pressure  and  locating  the  point  at  which  the  curve  is 
double  valued.  That  is,  there  are  two  densities,  corresponding  a  rarefied  phase 
and  a  dense  phase,  that  have  the  same  pressure  and  minimal  free  energy.  The 
critical  temperature,  he,  can  be  found  by  finding  the  isothermal  pressure  curve 
that  has  an  inflection  point 

dp{d,  he) 


dd 


=  0 


(96) 
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or 


(97) 


^  2rd{l  —  d){l  —  2d) 

To  verify  our  theory,  we  can  perform  exact  numerical  simulations  of  a  finite 
temperature  multiphase  system.  We  can  extend  an  Appert-type  nonthermal 
model  to  work  in  a  finite  temperature  domain  by  coupling  the  long  range  in¬ 
teractions  to  a  heat  bath  of  variable  density,  denoted  by  a  parameter  h,  and 
by  allowing  repulsive  long  range  interactions  in  addition  to  the  attractive  ones. 
This  is  done  in  such  a  way  that  the  likelihood  of  an  attractive  and  repulsive 
interaction  goes  as  (1  —  and  /i^  respectively  [15].  Figure  1  depicts  the 
long-range  interaction  called  bounce-back  and  how  it  is  coupled  to  a  heat  bath. 
At  zero  temperature,  when  /i  =  0,  we  recover  the  minimal  model  as  only  at¬ 
tractive  long  range  interactions  can  occur.  Fiqure  2  shows  the  time  evolution 
of  the  phase  separation  process  in  this  case  at  a  density  d  =  0.07  and  inter¬ 
action  range  r  =  61.  As  h  increases,  the  likelihood  of  repulsive  interactions 
also  increases  to  the  point  where  at  the  infinite  temperature  limit,  h  =  ^,  the 
likelihood  of  attractive  and  repulsive  interactions  becomes  equal.  The  occur¬ 
rence  of  both  long  range  attractive  and  repulsive  interactions  is  identical  to  a 
system  with  finite-impact  parameter  collisions.  Therefore,  the  infinite  temper¬ 
ature  system  behaves  as  an  ideal  neutral  fluid  but  with  an  enhanced  mean-free 
path.  The  nominal  strength  of  the  interaction  when  coupled  to  a  heat  bath  is 
a{h)  =  — ao(l  —  +  cxah"^  =  — ao(l  —  2h)  given  a  local  momentum  change 

of  magnitude  6p  =  oiomc  due  to  long  range  interactions  of  range  r.  For  a  two 
dimensional  example  D  =  2,  we  use  the  following  values  m  =  c=  l  =  l,  B  =  6, 
ao  =  2.  With  the  fluid  at  rest,  the  pressure  is  then 


p  =  3d  —  3rd^{l  —  d)^(l  —  2h) 


and  the  critical  value  of  h  is 


he 


I 

2rd{l  —  d){l  —  2d) 


(98) 


(99) 


Figure  3  shows  liquid-gas  coexistence  curves  for  this  lattice-gas  model  at  three 
different  interaction  ranges:  r  =  7,  9,  and  11/.  Both  the  mean  field  theory 
calculation  and  the  exact  numerical  data  are  presented.  The  comparison  of  the 
theory  to  the  numerical  simulation  is  in  good  agreement.  In  the  Boltzmann  limit, 
the  probability  of  a  long  range  interaction  goes  as  =  (i^(l  — d)^.  It  is  expected 
that  this  estimate  which  neglects  all  particle  correlations  would  suffer  the  most 
at  low  densities  where  the  mean  free  path  between  local  collisions  becomes 
comparable  to  the  range  of  the  non-local  interactions.  This  may  account  for 
the  deviations  that  are  observed  at  low  densities.  The  mean  field  predictions  of 
the  critical  point  is  also  in  quite  good  agreement  with  the  numerically  obtained 
values.  Figure  4  shows  the  mean  field  calculation  and  exact  numerical  data 
taken  at  five  different  interaction  ranges:  r  =  7  through  11/.  The  calculated 
value  of  he  is  slightly  higher  then  the  measured  value  for  all  cases  indicating  a 
systematic  deviation. 
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For  the  liquid-gas  system,  the  dispersion  relation  (85)  reduces  to 


^  =  ±k^  1  -I-  - do{l  -  do){l  -  2do)Da{h)  ^  {vak)  ?  J D{r^k)-i'^^- . 

(100) 

Figure  5  shows  the  real  and  imaginary  parts  of  the  liquid-gas  dispersion  curve  for 
a  two  dimensional  system  with  a  density  of  d  =  0.20  and  momentum  exchanges 
of  Sp  =  —2mc  over  a  range  of  r  =  91.  Also  shown  for  comparison  purposes 
are  the  dispersion  curves  for  an  ideal,  viscous  fluid.  In  a  long-range  lattice-gas, 
since  the  kinematic  shear  viscosity  is  dependent  on  the  square  of  the  mean- 
free  path  which  in  turn  is  proportional  to  the  interaction  range,  the  following 
approximation  is  made  for  t'o 

Vo{r)  =  i^r=o  (^  +  0-^r'^)  ■  (101) 

Numerical  confirmation  of  the  parabolic  dependence  of  the  kinematic  shear  vis¬ 
cosity  on  the  interaction  range  is  presented  in  figure  6.  Several  lattice-gas  models 
were  tested  by  varying  the  strength  and  number  of  interactions.  Viscosity  mea¬ 
surement  were  made  by  the  method  of  a  decaying  sinusoid  and  were  done  for 
systems  at  a  density  of  60%  filling. 

10  Crystallization 

Introduced  in  this  section  is  a  lattice-gas  automaton  with  multiple  fixed-range 
interactions  that  possesses  a  liquid-solid  phase  transition.  In  the  previous  sec¬ 
tion,  we  have  tested  our  formalism  that  models  interparticle  potentials  in  the 
coarse-grain  limit  by  using  a  single  anisotropic  fixed-range  interaction  in  the 
lattice-gas  dynamics  for  discrete  momentum  exchange  between  particles  in  the 
microscopic  limit.  Here  a  direct  generalization  to  the  finite  temperature  liquid- 
gas  model  is  introduced  using  long-range  repulsive  and  attractive  interactions 
over  multiple  ranges.  For  crystallization  to  occur,  at  least  two  interaction  ranges 
are  necessary:  an  attractive  short-range  interaction  and  a  longer-range  repulsive 
interaction  resulting  in  a  kind  of  Wigner  crystal. 

10.1  New  Way  for  Molecular  Dynamics  Modeling 

To  model  a  more  realistic  crystal,  that  is  one  that  can  undergo  rigid-body  motion 
such  as  rotation  and  that  can  have  well  defined  edges  or  surfaces,  more  then 
two  interaction  ranges  are  required.  Usually  four  to  eight  interaction  ranges  are 
used  to  produce  a  Leonard- Jones  type  molecular  potential. 

The  shortest-range  interaction  creates  a  potential  well  that  stably  traps  a 
group  of  lattice-gas  particles.  This  group  of  particles  remains  in  a  localized 
configuration  and  behaves  as  a  single  collective  entity.  This  persistent  collective 
entity  is  referred  to  here  as  an  “atom” .  As  in  the  liquid-gas  system,  each  lattice- 
gas  particle  possesses  a  discrete  field  that  acts  along  the  lattice  directions.  But 
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now  as  many  lattice-gas  particles  are  grouped  together,  in  the  coarse-grained 
limit  they  act  as  a  single  particle  with  a  continuous  field  around  it.  It  can  behave 
like  a  charged  particle  and  repel  other  such  atoms  in  the  system  or  can  behave 
like  a  Leonard- Jones  particle  and  attract  other  atoms  depending  on  the  chosen 
interactions.  Starting  from  a  uniformly  random  configuration  at  d  =  0.1,  the 
lattice-gas  spontaneously  crystallizes  into  arrays  of  these  atoms.  The  emergent 
crystalline  lattice  is  hexagonal-close-packed.  A  two  dimensional  example,  with 
an  underlying  512  x  512  lattice,  of  this  time-dependent  crystallization  process 
is  given  in  figure  7.  The  resulting  crystal  is  in  a  hexagonal-close-pack  configura¬ 
tion  since  we  have  strived  to  make  the  coarse-grained  interatomic  potential  be 
radially  symmetric.^  Three  dimensional  512^  simulations  of  the  crystallization 
were  also  carried  out,  see  figure  10. 

It  is  possible  to  measure  the  density  cross-section  for  the  crystal  in  its  final 
equilibrium  state,  see  figure  8.  With  a  principle  crystal  direction  aligned  parallel 
with  the  x-axis,  average  density  cross-sectional  data  was  taken  for  a  512  x  512 
system;  that  is,  512  samples  were  averaged.  In  this  case,  the  lattice-gas  model 
had  six  interaction  ranges:  r  =  —2  ,— 7,  19,  21,  —24,  —26.  Here  the  negative 
sign  preceding  the  range  denotes  an  attractive  interaction  at  that  range.  The 
averaged  cross-section  data  very  closely  produces  a  Gaussian  shaped  curve. 

In  the  two  dimensional  numerical  simulation,  to  obtain  isotropy  in  the  macro¬ 
scopic  limit,  12  directions  are  used  for  long-range  momentum  exchanges  instead 
of  6.  This  is  possible  because  the  underlying  triangular  lattice  has  6  momen¬ 
tum  states  and  the  total  possible  number  of  central-body  momentum  exchange 
directions  is  always  twice  the  lattice  coordination  number.  With  12  momentum 
exchange  directions,  the  crystal  is  stable  under  translation  along  any  direction 
and  in  fact  can  undergo  free  rotation.  Therefore,  the  crystal  acts  very  much  like 
a  solid  rigid  body.  This  rigid  body  can  also  support  elastic  waves — shear  waves 
and  compressional  waves  have  been  observed. 

The  local  stability  analysis  of  the  equations  of  motion  for  the  system’s  linear 
response  as  as  carried  out  in  section  §8  is  directly  applied  to  this  case.  In  the 
short-wavelength  limit,  the  dispersion  is  identical  to  that  for  an  ideal,  viscous 
fluid.  However,  for  small  wavelengths,  there  is  a  crucial  difference,  see  figure  9. 
The  imaginary  part  of  the  dispersion  relations  has  a  positive  peak  at  about 
k  =  0.08.  This  implies  an  instability  in  the  lattice-gas  system  that  ultimately 
gives  rise  to  the  crystalline  structure  characterized  by  cell  size  Therefore, 
the  linear  response  calculation  gives  a  nearly  quantitative  prediction  about  the 
size  of  the  emergent  crystal’s  cell  size.  The  interaction  ranges  used  in  the  lin¬ 
ear  response  calculation  are  r  =  7,  19,  21,  26  with  corresponding  interaction 
strengths  a  =  —2,  2,  2,  —2,  with  density  d  =  0.1.  The  dashed  curves  are  the 
dispersion  relations  for  an  ideal,  incompressible,  viscous  fluid  presented  here  for 
comparison  purposes. 

^If  the  density  of  the  system  is  increased,  one  does  observe  a  transition  from  a  hexagonally- 
ordered  bubble  phase  to  ordered  and  random  stripe  phases.  In  the  context  of  lattice-gases, 
Rothman  has  shown  some  pictures  similar  to  figure  7  in  an  two-component  immiscible  lattice- 
gas  with  a  short  range  attractive  interaction  and  a  longer  range  repulsive  interaction  [11]. 
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10.2  The  Crystal  Reconfiguration  Process 

An  expected  phenomenon  that  occurs  in  the  early  stages  of  the  crystal  formation 
is  the  emergence  of  grain  boundaries  and  defects.  Over  time,  given  the  inherent 
fluctuations  of  the  lattice-gas  dynamics,  the  crystal  undergoes  an  annealing 
process  that  removes  the  defects  and  eventually  produces  a  prefect  crystal.  In 
the  two-dimensional  case  with  a  radially  symmetric  coarse-grained  potential, 
the  hexagonal-close-pack  crystal  structure  emerges  as  just  mentioned.  Defect 
pairs  with  five  and  seven  neighbors  are  observed. 

An  unexpected  phenomenon  that  occurs  in  the  exact  simulation  is  the  process 
by  which  a  defect  is  removed.  To  describe  this  process,  consider  for  example,  an 
atom  with  five  neighbors.  It  may  persist  in  such  a  frustrated  situation  for  some 
time.  Yet  what  eventually  occurs  is  that  the  lattice  structure  near  this  defect 
begins  to  fluctuate — “tremors”  in  the  crystal  structure  are  observed.  That  is, 
the  other  atoms  in  the  immediate  vicinity  of  the  defect  begin  to  vibrate  about 
their  metastable  positions.  The  magnitude  of  the  fluctuations  increases  over 
time.  In  fact,  the  magnitude  of  the  fluctuations  appears  to  grow  and  one  may 
even  say  that  the  temperature  of  the  local  crystalline  structure  appears  to  rise. 
When  a  high  enough  local  temperature  of  this  sort  is  reached,  the  microscopic 
dynamics  suddenly  reconfigures  a  cluster  of  the  atoms  and  the  defect  vanishes. 
The  reconfiguration  of  the  atoms  usually  entails  a  small  local  rotation  of  a  cluster 
of  atoms  surrounding  the  defect  site. 

One  way  to  characterize  the  fluctuations  that  occur  in  the  system  leading  to 
the  rather  sudden  reconfiguration  of  a  cluster  of  atoms  is  to  compare  the  state 
of  the  system  at  one  time,  t,  to  the  state  at  some  later  time,  t  +  T,  by  computing 
a  hamming  length.  In  discrete  models,  such  as  the  Ising  model  or  a  Hopfleld 
neutral  network,  a  hamming  length  is  well  defined.  In  a  lattice-gas,  one  can  also 
define  a  hamming  length,  hi,  and  in  particular  we  do  so  in  the  coarse-grained 
limit.  That  is,  a  block  average  of  the  lattice-gas  number  variables  is  taken  to 
determine  a  density  held,  p{x).  The  hamming  length  is  calculated  by  summing 
over  all  points  of  the  density  held  as  follows 

hi  =  ^e{\p{x,t  +  T)  -  p{x,t)\  -  e) .  (102) 

X 

where  0{y)  is  the  step  function  which  is  zero  for  negative  y  and  unity  for  positive 
y  and  where  e  is  a  small  threshold  value.  Fiqure  11  shows  a  time  series  of  the 
hamming  length  for  a  long-range  lattice-gas  of  the  type  described  above.  The 
lattice  size  was  1024  x  1024,  block  size  used  was  8x8,  and  the  sampling  time 
was  T  =  10.  The  reason  for  measuring  the  hamming  length  is  that  it  provides 
a  rather  direct  and  simple  way  of  determining  the  scale  of  the  domain  of  atoms 
that  participate  in  a  crystal  reconfigure  event.  It  is  interesting  to  And  that 
the  domain  sizes  of  these  reconfigures  shows  power  law  behavior,  see  figure  12 
where  give  a  log-log  plot  of  the  frequency  of  occurance  of  a  reconfiguration 
versue  its  hamming  length.  In  figure  12  we  see  a  peak  at  hi  ~  5,  which  is  the 
most  common  background  fluctuation.  Larger  scale  fluctuations  occur  but  the 
probability  of  occurance,  p,  clearly  drops  off  according  to  a  power  law  of  the 


19 


form  p  oc  t4^.  In  this  case  alpha  is  approximately  6.3.  Smaller  fluctuations 

also  occur  but  these  are  not  responsible  for  the  reconfiguration  events  observed 
during  crystallization. 


11  Conclusion 

We  have  presented  a  mean-field  theory  of  lattice-gases  with  long-range  interac¬ 
tions.  We  have  focused  on  central  body  interactions  that  are  mediated  by  mo¬ 
mentum  exchange  events  between  remote  spatial  sites  and  have  used  these  type 
of  interactions  to  model  two  types  of  physical  systems:  a)  a  finite-temperature 
liquid-gas  dynamical  system;  and  b)  a  solid-state  molecular  dynamical  system. 
The  latter  lattice-gas  model  is  very  compelling  and  is  the  most  important  result 
of  this  paper.  This  lattice-gas  model  of  a  crystallographic  solid-body  offers  an  al¬ 
ternative  to  traditional  molecular  dynamics  modeling.  The  dynamical  behavior 
of  the  lattice-gas  solid  is  exactly  computed,  in  that  there  is  exact  conservation 
of  mass,  momentum,  and  energy.  A  solid  phase  is  self-consistently  produced 
through  the  collective  and  non-linear  behavior  of  billions  of  lattice-gas  particles 
as  they  interact  via  local  collisions  and  long-range  interactions.  A  linear  stability 
analysis  is  presented  that  predicts  the  formation  of  “molecules”  of  a  character¬ 
istic  intermolecular  spacing,  each  molecule  itself  occupying  a  finite  volume  and 
composed  of  thousands  of  lattice-gas  particles.  Therefore,  the  molecule,  is  not 
a  point  particle  but  is  distributed  over  several  lattice  sites  and  is  stable  in  a 
self-consistent  way.  Each  molecule  possesses  a  Lenord-Jones  type  potential  in 
the  coarse-grained  limit.  The  mass  of  a  molecule  as  well  as  its  held  are  both 
manifestations  of  the  spatial  distribution  of  lattice-gas  particles.  A  lattice-gas 
particle  is  interpreted  as  an  informational  token  that  composes  only  a  small  piece 
of  the  molecule  and  contributes  to  a  small  piece  of  its  held.  We  have  observed 
an  annealing  process  where  defects  are  removed  from  the  crystal  where  there  is 
a  succession  of  localized  vibrations  that  continually  build  to  the  point  where  a 
cluster  of  molecules  around  the  defect  can  eventually  undergo  a  reconfiguration. 
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Lattice  Multipole  Theorem 
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Theorem  1  Let  F  be  a  scalar  function  and  r  he  a  fixed  scalar  distance.  Then 
^aiF{x  ±  rea)  may  he  expressed  as  a  perfect  gradient  of  a  series  expansion: 


eaiF{x  ±  rea)  =  ±di 

a 


lD{rd)F{x)  , 


(103) 


where  Iu{z)  is  the  hyperbolic  Bessel  function. 


Proof:  We  begin  by  expanding  F{x  ±  rba)  as  follows 

^eaiF{x±re^)  =  (104) 

a  a 

=  eaz  [cosh.{rcreajdj)  ±  sinh(ro.eai9j)]  F{x). 

a 

(105) 

Now  the  directional  sum  over  Cai  cosh^ra-eajdj)  must  vanish  since  this  is  an  odd 
function  in  the  lattice  vectors.  Therefore,  the  expression  for  F(x  Frea)  reduces 
to 


Y  ±  rea) 

=  ±  X!  Ga*  s\vF{r^eajdj)F{x) 

(106) 

OO  2n— 1  2n— lQ2n— 1 

-  ±Ee-E  (2n-l)! 

a  n—l  ^  ' 

(107) 

°°  „2n-l  /  \ 

(108) 


Using  the  property  of  the  isotropic  lattice  tensors  of  rank  2n,  identity  (8),  and 
the  fact  that  the  total  number  of  terms  comprising  the  tensor  is  (2n—  1)!!, 

we  can  evaluate  the  directional  sum 


eaiF{x  ±  rea) 


n—l 


^2n-l 


(2n  -  1)!  D{D  +  2)  •  •  •  (D  +  2n  -  2) 


OO 


(2n 


2)!!£i(L»  +  2)---(£»  +  2n-2)' 


(109) 

(110) 


This  gives  F{x±  re^)  as  a  perfect  gradient 


^  ^  ^aiF(^X  dz  CGa)  —  FOi 


rB  Y: 


n—0.. 


(2n)!!(2n  +  D)!! 


F{x) 


or,  explicitly  writing  out  the  lowest  order  terms. 


Y  eaiF{x  ±  rea)  =  ±di 


rB 


r-^B 


D  ^  2D{D  +  2) 
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(111) 


(112) 
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Lemma  1  Let  D  he  a  positive  integer  and 

(113) 

so  that  e  is  zero  for  even  D  and  unity  for  odd  D.  Then  the  following  identity 
for  the  gamma  function  holds 


r 


+  1  + 


{2n  +  D)ll_. 

n  L  -  ^ 

9"+-^ 


(114) 


Proof  of  Lemma  1;  This  identity  for  the  gamma  function  follows  directly  by 
manipulating  r(n  +  1  +  ^)  separately  for  the  cases  when  D  is  even  and  odd. 
The  following  two  identities  for  the  gamma  function  are  useful 


r(n+l) 


nr(n)  =  n\ 

pOO 

2  /  dt  =  \/tt, 

Jo 


(115) 

(116) 


the  first  identity  holding  true  for  integer  n.  Using  these  identities  iteratively 
the  following  expression  holds 


r 


(2n-l)!!  /n 

2"  1^2  J 


(117) 


Then  generalizing  this  last  expression,  the  proof  of  the  lemma  follows  directly. 
QED 

Using  Lemma  1,  we  can  substitute  the  gamma  function  in  place  of  the  (2n  + 
D)ll  term  in  the  denominator  of  (111),  and  using  (2n)!!  =  2"n!,  we  then  have 
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Now  the  hyperbolic  Bessel  function  has  an  identical  series  expansion 


{-br 

n!r(n  +  \  +  v) 


(119) 


Therefore,  taking  v  =  ^  and  considering  z  to  be  a  differential  operator,  z  — >■  rd, 
then  (118)  becomes 


eaiFix  ±  roa)  =  ±rdi 


’iloird)Fix) 


(120) 
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which  completes  the  proof  of  Theorem  1.  QED 

Numerical  Implementation  of  the  Long-Range  2-Body  Interactions 

A  simple  computational  scheme  is  employed  that  allows  all  the  dynamics  to 
be  computed  in  parallel  with  two  additional  bits  of  local  site  data,  for  outgoing 
and  incoming  messengers,  regardless  of  the  number  of  long-range  neighbors. 
The  computational  scheme  is  an  efficient  decomposition  of  a  lattice-gas  with 
many  neighbors.  It  is  conceptually  similar  to  the  idea  of  virtual  intermediate 
particle  momentum  exchanges  that  is  well  known  in  particle  physics.  All  2- 
body  interactions  along  a  particular  direction  define  a  spatial  partition  that  is 
updated  in  parallel.  Random  permutation  through  the  partitions  is  sufficient  to 
recover  the  necessary  isotropy  as  long  as  enough  momentum  exchange  directions 
are  used. 

An  interparticle  potential,  y(x  —  x'),  acts  on  particles  spatially  separated 
by  a  fixed  distance,  x  —  x'  =  r.  An  effective  interparticle  force  is  caused  by  a 
non-local  exchange  of  momentum.  Momentum  conservation  is  violated  locally, 
yet  it  is  exactly  conserved  in  the  global  dynamics. 

For  the  case  of  an  attractive  interaction,  there  exists  a  bound  states  in  which 
two  particles  orbit  one  another.  Since  the  particle  dynamics  are  constrained  by 
a  crystallographic  lattice  we  expect  polygonal  orbits.  In  figure  13  we  have  de¬ 
picted  two  such  orbits  for  a  hexagonal  lattice-gas.  The  range  of  the  interaction  is 
r.  Two-body  single  range  attractive  interactions  are  depicted  in  figures  13b  and 
13c,  the  bounce-back  and  clockwise  orbits  respectively.  Momentum  exchanges 
occur  along  the  principle  directions.  The  interaction  potential  is  not  spherically 
symmetric,  but  has  an  angular  anisotropy.  In  general,  it  acts  only  on  a  finite 
number  of  points  on  a  shell  of  radius  | .  The  number  of  lattice  partitions  neces¬ 
sary  per  site  is  half  the  lattice  coordination  number,  since  two  particles  lie  on  a 
line.  Though  microscopically  the  potential  is  anisotropic,  in  the  continuum  limit 
obtained  after  coarse-grain  averaging,  numerical  simulation  indicates  isotropy 
is  recovered. 


A  Simple  Example:  Bounce-Back  Orbit 

A  long-range  lattice-gas  of  the  type  we  are  considering  still  possesses  the  usual 
local  dynamics  of  a  hydrodynamic  lattice-gas.  To  extend  the  local  lattice-gas 
update  rules  to  include  long-range  interactions,  we  use  two  additional  bits  of 
local  site  data.  This  will  allow  us  to  implement  a  long-range  interaction  using 
strictly  local  updating  and  therefore  our  algorithm  will  still  be  parallelizable 
just  as  a  usual  local  lattice-gas.  The  two  additional  bits  will  denote  the  occupa¬ 
tion  numbers  of  messenger  particles,  or  “photons” .  The  idea  of  using  messenger 
particles  was  introduced  by  Appert  et  al.[l].  We  have  two  types  of  messen¬ 
ger  states,  to  represent  incoming  and  outgoing  conditions,  and  we  denote  the 
messengers  as  ijji  and  ipr- 

Now  for  the  simplest  long-range  lattice-gas  model,  we  therefore  use  eight  bits 
of  local  site  data.  Since  long-range  interactions  occur  between  remote  spatial 
sites,  say  x  and  x' ,  the  messenger  particles  will  travel  either  parallel  or  antipar- 
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allel  to  the  vector  r  =  x  —  x'.  All  pairs  of  sites  throughout  the  entire  space  that 
are  separated  by  vector  r  can  therefore  all  be  updated  in  parallel.  We  refer  an 
update  step  of  all  pairs  of  2-body  interactions  along  direction  r  as  a  partition, 
denoted  by  r^..  All  possible  two-body  interaction  pairs  are  then  computed  by 
performing  all  possible  partitions  of  the  space.  So  it  requires  many  scans  for 
the  space  to  perform  a  single  long-range  interaction  step. 

In  our  two-dimensional  example  using  a  triangular  lattice,  there  are  three 
possible  partitions.  The  number  of  partitions  is  never  smaller  than  half  the 
lattice  coordination  number.  In  the  two-dimensional  case,  the  simplest  long- 
range  lattice-gas  algorithm,  though  perhaps  not  the  most  efficient  algorithm,  is 
to  use  three  sequential  scans  of  the  space,  each  scan  performing  the  updating 
necessary  for  a  single  partition,  see  figure  13a.  Often  times,  depending  on 
the  complexity  of  the  long-range  interactions  and  the  dimensionality  of  the 
lattice,  it  is  possible  to  perform  simultaneous  updating  of  multiple  partitions. 
This  of  course  is  more  desirable  yet  causes  more  complexity.  Furthermore,  this 
updating  requires  an  extra  pair  of  messenger  particles  for  each  partition  to  be 
simultaneously  updated.  For  simplicity,  we  will  not  deal  with  this  case  here, 
however  our  implementation  on  the  CAM-8  does  use  simultaneous  partition 
updating — repulsive  and  attractive  partitions  are  performed  in  parallel  using 
four  messenger  bits. 

Let  us  consider  a  simple  example  of  the  long-range  lattice-gas  algorithm, 
the  minimal  model  of  Appert.  Here  we  consider  only  bounce-hack  attractive 
interactions.  Suppose  there  is  a  single  particle  at  site  x  =  0  and  there  is  also 
a  single  particle  at  site  x'  =  ri;  that  is,  no(x)  =  1,  n3(x)  =  0,  no(x')  =  0  and 
n3(x')  =  1  with  all  other  bits  at  x  and  x'  being  zero,  see  figure  13b.  Here  we  are 
using  the  bit  convention  shown  in  table  1.  Then  the  two  particles  are  separated 
by  a  distance  r  and  are  moving  away  from  each  other.  The  attractive  long-range 
interaction  will  effectively  flip  their  respective  directions  making  no(x)  =  0, 
n3(x)  =  1,  no(x')  =  1  and  n3{x')  =  0  so  that  the  two  particles  will  now  be 
moving  toward  each  other.  There  is  a  local  momentum  change  of  2mci  at  x' 
and  an  opposite  momentum  change  of  — 2mci  at  x.  Locally  momentum  is  not 
conserved,  but  nonlocally  it  is. 

The  first  step  of  the  long-range  interaction  sequence  is  to  choose  a  partition, 
say  Fr,  and  then  to  emit  messenger  particles  along  the  partition  axis.  The  basic 
local  rule  for  this  first  step  is  the  following:  a  photon  is  emitted  at  a  site  if 
there  exists  a  particle  at  that  site  that  can  partake  in  a  long-range  interaction. 
Another  way  of  expressing  this  rule  is:  send  only  if  you  can  receive.  Obviously, 
for  a  particle  to  partake  in  an  interaction  there  must  be  both  a  particle  and  a 
hole  at  that  site.  The  factorized  probability  of  having  such  a  situation  is  just 
d{l  —  d).  So  to  continue  with  our  example,  for  a  photon  to  be  emitted  at  some 
site  X  parallel  or  antiparallel  to  a  partition  direction  i,  we  use  the  following  rule 

xfj.{x)  =  no(x)(l  -  n3(x))  (121) 

iii{x)  =  n3(x)(l  -  no(x)).  (122) 

Note  that  according  to  this  local  rule,  only  one  photon  can  be  created  at  a  site. 
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and  consequently  we  eliminate  the  possibility  of  a  long-range  interaction,  say  of 
range  2r,  mediated  through  a  doubly  occupied  site.  The  important  consequence 
of  the  emission  step,  is  that  for  two  sites  separated  by  the  interaction  distance, 
r,  if  both  sites  send  photons,  both  will  necessarily  receive  them,  which  strictly 
enforces  nonlocal  momentum  conservation.  Give  and  ye  shall  receive  (provided 
your’s  is  received).  Letting  ipa  =  ipr  and  pJ-a  =  ipi,  in  general  we  can  write  the 
emission  step  of  the  minimal  interaction  as 

ipa{x)  =  n_a{x){l  -  na{x)),  (123) 

where  a  =  0, 1,  2  covers  all  the  partitions. 

After  the  emission  step,  follows  a  long-range  kick  of  the  messenger  bits.  In 
the  simple  example,  all  photons  ipi  are  kicked  along  — ri  and  all  photons  ipr  are 
kicked  along  ri.  In  general  for  the  long-range  kick  we  have 

tp'^  {x  +  rea)=  pJa  (^)  ■  (124) 

Finally,  we  have  the  absorption  step  of  the  long-range  interaction  sequence. 
Here  the  local  particle  momentum  state  is  updated  as  the  particles  flip  their 
directions  in  our  example 

n^{x)  =  nsix)  + 'tp'i{x)no{x){l  -  nsix))  - 'ipP{x)n3{x){l  -  no{x)) 

(125) 

no(f)  =  no{x)  +  'ipP{x)n3{x){l-no{x))-'ip'i{x)no{x){l-n3{x)). 

(126) 

Moreover,  in  this  step  all  the  messenger  bits  are  set  to  zero  throughout  the 
entire  space.  For  any  direction,  the  local  absorption  rule  could  more  simply  be 
written  as 

n'aix)  =  na{x)  +  'tp'_^{x)tpa{x)  -  tp'^{x)tp  -  a{x)  .  (127) 

Substituting  in  (123)  and  (124)  into  (127),  we  have  a  single  boolean  expression 
in  terms  of  number  variables  for  a  single  long-range  interaction  step  for  partition 
Fr  as  follows 

n^{x)  =  na{x)  + 

na(x  +  reo)(l  -  n_a(x  -I-  rea))n_a(x)(l  -  na{x))  - 
n-a(x  -  rea)(l  -  na(x  -  rea))na(x)(l  -  n_a(f)) 

(128) 

For  convenience  we  define  a  long-range  collision  operator.  Pa,  as  follows 

Pa{x,  X  +  rCa)  =  1p-a{x  +  rea)tpa{x)  =  1p'_a{x)lpa{x) ,  (129) 

so  that  we  may  write 

n'aix)  =  na{x)  +  Pa{x,x  +  Tta)  -  P-a{x ,  X  -  rCa)  ■  (130) 

The  state  data  for  this  simple  example  we  have  been  considering  is  given  in 
table  2,  which  represents  all  the  steps  of  a  long-range  interaction  sequence  for 
a  partition  along  the  x-axis. 
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B  Another  Example:  Clockwise  Orbit 

To  continue  illustrating  our  implementation  of  a  long-range  lattice-gas,  in  this 
section  we  again  consider  a  system  with  a  single  attractive  interaction  of  range 
r,  however  the  local  momentum  states  participating  in  the  interaction  are  not 
along  the  partition  direction.  Yet  in  the  example  given  here,  the  momentum 
exchange  is  still  along  the  partition  direction  so  that  the  long-range  interac¬ 
tion  remains  a  central-body  one,  resulting  in  a  bound  state  with  two  particles 
trapped  in  a  clockwise  orbit.  (Note  that  the  restriction  to  central-body  forces 
is  not  necessary,  but  is  presented  here  for  convenience.)  In  this  slightly  more 
complicated  example,  the  local  rules  for  photon  emission,  and  absorption,  (123) 
and  (127)  respectively,  have  a  more  general  form  with  the  implication  that  the 
emission  and  absorption  of  photons  is  different  from  the  previous  example  of 
the  bounce-back  orbit  and  should  be  noted  when  making  look-up  tables  to  do 
this  computation.  The  local  photon  emission  rules  can  be  written 

tpaix)  =  nc{x){l  -  nd{x))  (131) 

ip-aix)  =  ng{x){l  -  nh{x))  (132) 

where  the  bits  c,d,g,h  must  by  chosen  so  momentum  is  conserved 

Sc  -  ed  +  eg  -  Ch  =  0  (133) 

as  well  as  be  constrained  by  central-body  parallel  and  perpendicular  momentum 
exchange  conditions 


{ec- ed- eg  +  Ch)  ■  r  =  2Ap  (134) 

{ec- ed- eg  +  Ch)  X  r  =  0,  (135) 


where  Ap  is  the  momentum  change  per  site  due  to  the  long-range  interaction. 
In  (131)  and  (132)  the  difference,  over  our  previous  example  of  the  bounce-back 
orbit,  is  the  possibility  of  having  two-photons  emitted  at  a  single  site. 


To  be  explicit,  for  the  two-dimensional  triangular  lattice,  we 
(133),  (134),  and  (135)  by  choosing  the  indices  c,d,g,h  as  follows: 

can  satisfy 

c  =  a  —  2 

(136) 

d  =  a  —  1 

(137) 

g  =  -c 

(138) 

h  =  —d. 

(139) 

An  example  of  this  choice  of  indices  is  illustrated  in  figure  13c.  Then  the  emission 
rule,  (131)  and  (132),  is  simply 

1pa{x)  =  na-2{x){l  -  na-l{x))  (140) 

Since  the  kicking  of  the  photons  is  the  same  in  this  example  as  in  the  previous 
one,  (124)  still  holds 

+  rta)  =  -tpa^x). 
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By  re-expressing  (127)  more  generally,  we  can  write  a  local  absorption  rule 


<(^)  =  na(x)  +  V'-(a-Hl)(^?)V'a-Hl(^)  “  V'l-l  (^)V'-(a-l)  (^)  (141) 

or  more  elegantly 

n(j(f)  =  na{x)  +  Pa+i{x,x  +  rea+i)  -  P-a+i{x,x  -  rca-i).  (142) 

Substituting  in  (140)  and  (124)  into  (141)  and  after  some  manipulation  of  the 
indices,  we  have  a  single  boolean  expression  in  terms  of  number  variables  for  a 
single  long-range  interaction  step  for  partition  as  follows 

n'aix)  =  na{x)  + 

na+2{x  +  rea+i)(l  -  n_a(f  -I-  rea+i))na-i(f)(l  -  na{x))  - 

n-a{x  -  rea-i){l  -  na-2{x  -  rea-i))na{x){l  -  na+i{x)). 

Table  3  gives  the  local  site  data  for  the  x-axis  partition  of  a  clock- wise  orbit. 
The  particle  n4{x)  acts  as  a  kind  of  spectator  is  this  example,  illustrating  that 
two  photons  can  be  emitted  from  a  single  site.  It  is  also  possible  to  have  two 
photons  absorbed  at  a  single  site.  Let  us  consider  a  back-to-back  interaction  over 
three  sites.  Suppose  there  are  particles  at  sites  x  =  0,  x'  =  ri,  and  x"  =  2ri. 
Table  4  gives  the  site  data  for  these  sites  were  there  are  two  photons  emitted 
and  absorbed  at  x'  in  the  middle  location. 
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Table  2:  Long-range  interaction  sequence 


events 

na{x) 

Zl{x) 

Zr{x) 

na{x') 

Zl{x') 

Zr{x') 

initial 

100000 

0 

0 

000100 

0 

0 

emit 

100000 

0 

1 

000100 

1 

0 

kick 

100000 

1 

0 

000100 

0 

1 

absorb 

000100 

0 

0 

100000 

0 

0 

Table  3:  Long-range  interaction  sequence  with  two  photons  emitted  at  a  single 
site 


labels 

naix) 

Zlix) 

Zr{x) 

na{x') 

Zl{x') 

Zr(x') 

events 

010010 

0 

0 

000010 

0 

0 

emit 

010010 

1 

1 

000010 

1 

0 

kick 

010010 

1 

0 

000010 

0 

1 

absorb 

001010 

0 

0 

000001 

0 

0 

Table  4:  Long-range  interaction  sequence  with  two  photons  emitted  and  ab¬ 
sorbed  at  site  x'  in  a  back-to-back  interaction 


events 

na(x) 

Zl{x) 

Zr(x) 

na(x') 

zi(x') 

Zr(x') 

na(x") 

zi(x") 

Zrix") 

initial 

010000 

0 

0 

010010 

0 

0 

000010 

0 

0 

emit 

010000 

0 

1 

010010 

1 

1 

000010 

1 

0 

kick 

010000 

1 

0 

010010 

1 

1 

000010 

0 

1 

absorb 

001000 

0 

0 

001001 

0 

0 

000001 

0 

0 
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Figure  1:  Long-range  bounce-back  collision  coupled  to  a  heat  bath.  Hallow  circle  denotes 
a  heat-bath  hole  and  a  filled  circle  denotes  a  heat-bath  particle.  Interaction  (a)  to  (b)  repre¬ 
sents  particle  attraction  emitting  two  heat  bath  particles.  The  reverse  interaction  (b)  to  (a) 
represents  particle  repulsion  absorbing  two  heat  bath  particles. 


Figure  2:  Time  evolution  of  a  liquid-gas  phase  separation  for  a  lattice  gas  with  long  range 
attractive  interactions  at  range  r  =  6  on  a  1024  x  1024  lattice  starting  with  a  uniformly 
random  configuration  of  density  d  =  0.07. 
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p  =  3  d  -  21  ci"'2  (  l-d)  -^2  (1-211) 


p  =  3  d  -  27  d-^2(l— d)-^2  (1— 21i) 


p  =  3  d  -  33  d-^2(l-d)''2  (l-21i) 


ir=  7 


jr=  9 


zr=ll 


Figure  3:  Liquid-gas  coexistence  curves  for  the  simplest  lattice  gas  models  (Appert-type 
minimal  model  extended  to  finite  temperatures  by  coupling  to  a  heat  bath  with  filling  fraction 
h)  with  a  long-range  attractive  and  repulsive  interactions.  Mean  field  theory  versus  numerical 
data  is  shown  for  the  model  at  three  different  interaction  ranges. 
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Figure  4:  Liquid-gas  critical  point  values  versus  interaction  range  for  the  simplest  lattice  gas 
models  (Appert-type  minimal  model  extended  to  finite  temperatures  by  coupling  to  a  heat 
bath  with  filling  fraction  h)  with  a  long-range  attractive  and  repulsive  interactions.  Mean 
field  theory  versus  numerical  data  is  shown  for  the  model  at  five  different  interaction  ranges. 
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Figure  5:  The  real  and  imaginary  parts  of  the  liquid-gas  dispersion  relation  for  a  2- 
dimensional  system  rendered  with  the  solid  curves.  Models  parameters:  range=9,  density=0.2, 
attractive  interaction  of  strength=2.  The  dashed  curves  are  the  dispersion  relations  for  an 
ideal,  incompressible,  viscous  fluid. 
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Kinematic  Viscosity 


Figure  6:  Kinematic  shear  viscosity  versus  nonlocal  interaction  ranges  for  several  liquid-gas 
lattice  gas  models  at  density  =  0.6. 


34 


t  =  50 


t  =  500 


t  =  5000 


Figure  7‘  Time  evolution  of  crystallization  in  a  2-dimensional  lattice  gas  with  multiple  fixed- 
range  2-body  interactions.  The  resulting  crystal  is  in  a  hexagonal-close-pack  configuration 
since  the  coarse-grained  interatomic  potential  is  radially  symmetric.  The  underlying  lattice  is 
512  X  512.  Started  with  a  uniformly  random  configuration  at  d  =  0.1.  Twelve  directions  are 
used  for  long-range  momentum  exchanges.  Grain  boundaries  and  defects  are  observed  during 
the  early  stages  of  the  crystal  formation. 
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Figure  8:  Average  density  cross-section  of  an  equilibrium  lattice  gas  crystal  formed  using 
six  interaction  ranges:  r  =  —2  7,  19,  21,  —24,  —26  (negative  ranges  denotes  attraction). 

A  string  of  Gaussian  functions  provides  an  excellent  fit  to  the  numerical  data  for  a  two 
dimensional  system. 
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Figure  9:  Rendered  the  solid  lines  are  the  real  and  imaginary  parts  of  the  crystal  dispersion 
relation  for  for  small  wavelengths  for  a  2-dimensional  system.  Models  parameters:  ranges=7, 
19,  21,  26,  density=0.1,  interaction  of  strength=-2,  2,  2,  -2.  The  dashed  curves  are  the 
dispersion  relations  for  an  ideal,  incompressible,  viscous  fluid.  A  positive  peak  is  observed 
in  the  imaginary  part  of  the  dispersion  relation  giving  rise  to  an  unstable  growth  of  small 
perturbations  causing  an  emergent  crystal  structure. 
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Figure  10:  Lattice  gas  three-dimensional  crystal.  The  underlying  lattice  is  a  512  cube.  The 
simulation  was  done  on  a  128-million  node  CAM-8. 
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Figure  11:  Hamming  length  time  series  data  for  a  long-range  lattice  gas  with  8  different 
ranges.  The  hamming  length  corresponds  to  the  size  of  a  reconfiguration  event  due  to  inherent 
fluctuation  of  the  underlying  lattice  gas.  Data  taken  for  a  1024  x  1024  simulation  using  8x8 
block  averaging  and  a  sample  time  T  =  10. 
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Figure  12:  Power  law  behavior  of  the  frequency  of  occurrence  versus  hamming  length  or 
the  domain  size  involved  in  a  reconfiguration  event.  Graph  determined  from  hamming  length 
time  series  data  presented  above. 


Figure  13:  Simple  bound-state  orbits  due  to  a  long-range  attractive  interaction  where  the 
dotted  path  indicates  the  particle’s  closed  trajectory:  (a)  partition  directions;  (b)  bounce-back 
orbit  with  |Ap|  =  2  and  zero  angular  momentum;  and  (c)  clockwise  orbit  with  |Ap|  =  1  and 
one  unit  of  angular  momentum.  Head  of  the  dashed  arrows  indicates  particles  entering  the 
sites  of  partition  ro  at  time  t.  Tail  of  the  black  arrows  indicates  particles  leaving  those  sites 
at  time  t  t.  The  counter-clockwise  orbit  is  not  shown. 
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